KRMr from picture to TS

Sven Gastauer

2023-03-06

Data included in the package

The KRMr package includes example shape data for sardines found on the NOAA Southwest Fisheries Science Center webpage (SWFSC). The data contains the shape descirption for the fish body data(pil_fb) and the swimbladder data(pil_sb).

The shape can also be plotted:

library(KRMr)
library(knitr)

data(pil_sb)
data(pil_fb)
par(mfrow=c(1,2))
KRMr::shplot(x_fb = pil_fb$x_fb, w_fb = pil_fb$w_fb,
            z_fbU = pil_fb$z_fbU, z_fbL = pil_fb$z_fbL,
            x_sb = pil_sb$x_sb, w_sb = pil_sb$w_sb,
            z_sbU = pil_sb$z_sbU, z_sbL = pil_sb$z_sbL)

More example data will be added eventually.

Model input parameters

The KRMr implementation of the Kirchhoff Ray Mode model takes several input parameters that will define the resulting scattering. A minimal example can be used to run example simulations for sardines, with shape information found on the NOAA Southwest Fisheries Science Center webpage (SWFSC):

krm(
  frequency = 120 * 1000,
  c.w = 1490,
  rho.w = 1030,
  theta = 90,
  cs = c(1570, 345),
  rhos = c(1070, 1.24),
  shape = list(pil_fb, pil_sb),
  modes = c("fluid", "soft"),
  L = 0.21,
  fb = 1
)
## 2023-03-06 20:53:57:2 Body parts detected with 2 densities and 2 internal soundspeeds
##   frequency  c.w rho.w theta    L scale                        f1
## 1    120000 1490  1030    90 0.21     1 -0.001188948-0.000358712i
##                         f2        sigma       TS
## 1 0.001346062+0.000360853i 0.0001571285 -76.0749

Orientation

° The KRM model holds validity for orientations between 65° and 115°.
Orientation Definition

Shape definition

Shapes in the KRM model are scaled by the length parameter (in meters).
Shape Definition

Shapes are integrated over a number of elements along the shape, typically from tail to head. Each of these elements is described by \(x\) the position alog the x axis, along the body from tail to head. The lower (\(z_L\)) and upper points (\(Z_U\)) at a given \(x_i\), defining the height of the element, \(w\) is the width of the element.

Shapes will be combined coherently within the model. The model output will contain the complex signal from each shape and an incoherent sum can be computed ad hoc. Alternatively a coherent summation of selected parts could be computed as well.

Taking picture for shape extraction

To get good shape information from pictures or XRays, it is recommended to take pictures form the side and from the top of each shape the should be included in the model.
Here we illustrate how to take pictures for shape extraction through a few examples.

Using ImageJ to extract shapes

If you don’t have ImageJ, get it from the ImageJ webpage.
- Open the picture you want to process. - If the fish is not lying straight in the image (i.e. there is an angle at the head to tail line), then the image has to be aligned using Image->Transform->Rotate. To find the rotation angle, draw a line from the head to the tail and press M. This will open the Results Dialog and contain information on the angle of the line.
- If a ruler or another size reference is visible in the picture, draw a line across the target of known length. Click Analyze -> Set Scale and fill out the known length. This will calibrate the pixel size in the opened image. - Start drawing a polygon shape around the fish body or whichever part of the image you want to include in the model. Once completed, the XY information for the body part can be exported File -> Save as.. -> XY coordinates.

** Using the ImageJ ROI manager**

  • It is recommended to use thee ImageJ ROI manager. Open the ROI manager by clicking: Analyze -> Tools -> ROI Manager. Now you can add the shape selection by clicking Add. It is recommended to give the shape a sensible name by clicking Rename.

The procedure should be repeated for all relevant body parts. Once completed, the collection of ROIS can be exported, by selecting all shapes in the ROI manager, click More -> Save.

Read ImageJ ROI zip file into R as list:

library(RImageJROI)
side = RImageJROI::read.ijzip(fn) #fn should be the filename and path to the ROI.zip file
shapes = lapply(side, FUN=function(x){
  tmp=data.frame(x$coords)
  tmp$Part=x$name
  return(tmp)})

** Without ImageJ ROI Manager**
The XY Corrdinates exported from IMageJ can be read into R with read.table(). A 3rd column with a shape descriptor name has to be added. THis should be something snesible like ‘Dorsal_body’, ‘Dorsale_bladder’ or ‘Lateral_body’,‘Lateral_bladder’ etc.

Converting IMageJ XY Coordiantes into KRM shape

The KRMr function Imagej2shp(shp,dorsal = c("Dorsal_body", "Dorsal_bladder"),lateral = c("Lateral_body", "Lateral_bladder"), body = 1, xy = c(1, 2), nam = 3, n = 0.01) can be used to translate ImageJ coordinates into KRM shapes.

  • shp dataframe containing the shape information, with at least 3 colums, contianing x,y and descriptor of shape
  • dorsal names of shapes from dorsal aspect
  • lateral names of shapes from lateral aspect
  • body position of the largest body in the shape names, defaults 1
  • xy position of the xy coordinates in the dataframe defaults c(1,2)
  • nam position of name of the shape, defaults 3
  • n Resolution of the ouput shape where 1/n elements will be returned

The resulting shape can directly be used as shape input into the krm() function.

Example Pollock

Pollock Top Pollock Side
For the XRays, the fish body side view and top shape as well as the swimbladder side view and top view shapes were extracted. They were save as XY coordinates and pasted into one csv file with a third column containing the name of the shape.

We can read the shape information, translate it into a KRM shape and then display it in 3D:

shp1 = KRMr::Imagej2shp(shp=read.csv("Pollock01.csv"))
KRMr::get_shp3d(shp1)

Now we are ready to run a KRM simulation:

frequency=c(18,38,70,120,200)*1e3 #frequencies Hz
theta = seq(65,90,115) #orientations

c.w = 1490 #Ambient water m/s
rho.w = 1026.8 #density ambient water kg/m3

c.fb=1570 #soundspeed body 
c.b=345 #soundspeed swimbladder
cs = c(c.fb, c.b)

rho.fb=1070 #density body kg/m3
rho.b=1.24 #density swimbladder kg/m3
rhos = c(rho.fb, rho.b)

L = 29.7/1e2 #Lengths for pollocks in m

krm1 = krm(frequency = frequency,
      c.w = c.w,
      rho.w = rho.w,
      theta = theta,
      cs =cs,
      rhos = rhos,
      shape = shp1,
      modes = c("fluid","soft"),
      L = L[1],
      fb = 1)
knitr::kable(krm1)
frequency c.w rho.w theta L scale f1 f2 sigma TS
18000 1490 1026.8 65 0.297 0.297 -0.0011486-0.0001648i 0.0035069+0.0003175i 0.0023632 -52.52992
38000 1490 1026.8 65 0.297 0.297 -0.0008918-0.0011292i -0.0014728-0.0023896i 0.0042395 -47.45380
70000 1490 1026.8 65 0.297 0.297 0.0003944+0.0014709i 0.0000685+0.0010554i 0.0025684 -51.80687
120000 1490 1026.8 65 0.297 0.297 -0.0017041-0.0002549i 0.0012680+0.0020808i 0.0018773 -54.52950
200000 1490 1026.8 65 0.297 0.297 0.0016757-0.0002071i 0.0003115+0.0037990i 0.0041050 -47.73368

Example: Herring

Herring

Challenge: Only one top view is available for many Sideviews. Only side view of the swimbladder is avaialble.

For herring shown here, each element of the swimbladder is condiered to be circular (diameter = width for each point along the x axis (the body length).
For each fish, a picture should be taken sideways and from the top. The swimbladder should be clearly visible in the picture.

Here we use these images as an example to show how to use one image providing the view from the top to multile sidevidews and how to use only the sideview of the swimbladder to generate a swimbladder where each element is assumed to be circular (note that each element can have a different radius).

Read the available data:

fn = "HER.zip" # Path to the ImageJ ROI.zip file
top = read.table("HER_body_top.txt") # Read the only avaialble top view file

library(RImageJROI)
side = RImageJROI::read.ijzip(fn)
shapes = lapply(side, FUN=function(x){ #load all side view coordinates into a list
  tmp=data.frame(x$coords)
  tmp$Part=x$name
  return(tmp)})

WHat shape side views information is available?

 par(mfrow=c(3,4), mar=c(1,1,1,1)) #have a quick look at the available shapes
 for(i in seq(1,length(side), by=2)){
   plot(shapes[[i]][, 1],-shapes[[i]][, 2], asp=1, type="l")
   points(shapes[[i+1]][, 1],-shapes[[i+1]][, 2], asp=1, type="l")
 }

Now we have to add width information to each Herring:

####Add widths

names(top)=c("x", "y")
top$Part = "Body Top"
topr=data.frame(x=-top[,1], y=top[,2])
topr$Part = "Body Top"
shr = function(x){x$Part = "SwB Top";return(x)}
tops = list(top,shr(shapes[[2]]),  top,shr(shapes[[4]]),  top,shr(shapes[[6]]),  topr,shr(shapes[[8]]),
     topr,shr(shapes[[10]]),topr,shr(shapes[[12]]),topr,shr(shapes[[14]]),top,shr(shapes[[16]]),
     top,shr(shapes[[18]]), top,shr(shapes[[20]]), topr,shr(shapes[[22]]))

To add the top shape to all side views, we have to scale it correctly and then we can use the Imagej2shp() function:

#scale top
for(i in seq(1,length(side), by=2)){
  k=(i+1)/2
  message(Sys.time(), ": Processing shape ", k)
  xmin0 = min(shapes[[i]]$x); xmax0 = max(shapes[[i]]$x)
  dx0 = xmax0-xmin0
  ymin0 = min(shapes[[i]]$y)
  shapes[[i]]$x=(shapes[[i]]$x-xmin0)/dx0
  shapes[[i]]$y=(shapes[[i]]$y-ymin0)/dx0
  shapes[[i+1]]$x = (tops[[i+1]]$x - xmin0) / dx0
  shapes[[i+1]]$y = (tops[[i+1]]$y - ymin0) / dx0

  tmin=min(tops[[i]]$x); tmax=max(tops[[i]]$x)
  ymin = min(tops[[i]]$y)
  dx= tmax-tmin
  ty=min(tops[[i+1]]$x); tmax=max(tops[[i+1]]$x)
  tops[[i]]$x = (tops[[i]]$x - tmin) / dx
  tops[[i]]$y = (tops[[i]]$y - ymin) / dx

  tops[[i+1]]$x = (tops[[i+1]]$x - xmin0) / dx0
  tops[[i+1]]$y = (tops[[i+1]]$y - ymin0) / dx0

  shp = rbind(shapes[[i]],shapes[[i+1]],tops[[i]], tops[[i+1]])

  bi=sprintf("%02d", k)
  #write.csv(shp,paste0("HER_",bi,"_shp.csv"))

  shape=Imagej2shp(shp=shp,dorsal = c("Body Top","SwB Top") ,
                   lateral = c(unique(shp$Part)[1],unique(shp$Part)[2]),
                   body=1, xy=c(1,2), nam=3, n=0.01)
    saveRDS(shape,paste0("HER_",bi,"_shape.rds"))
}

Once completed we can create a subplot of all generated shapes in 3D:

### Generate 3D plots

library(plotly)
library(dplyr)

p3 = lapply(1:11,FUN=function(x) {
  message(paste0("scene",x))
  get_shp3d(readRDS(paste0("HER_",sprintf("%02d", x),"_shape.rds")), scene=paste0("scene",x))
  })
fig <- subplot(p3)%>% layout(scene=list(domain=list(x=c(0,0.4),y=c(0.75,1)),
                                        aspectmode='data'))

fig <- fig %>% layout(title = "Herring Shapes",
                      scene1 = list(domain=list(x=c(0,0.4),y=c(0.75,1)),
                                   aspectmode='data'),
                      scene2 = list(domain=list(x=c(0.3,0.7),y=c(0.75,1)),
                                    aspectmode='data'),
                      scene3 = list(domain=list(x=c(0.6,1),y=c(0.75,1)),
                                    aspectmode='data'),
                      scene4 = list(domain=list(x=c(0,0.4),y=c(0.5,.75)),
                                    aspectmode='data'),
                      scene5 = list(domain=list(x=c(0.3,0.7),y=c(0.5,0.75)),
                                    aspectmode='data'),
                      scene6 = list(domain=list(x=c(0.6,1),y=c(0.5,0.75)),
                                    aspectmode='data'),
                      scene7 = list(domain=list(x=c(0,0.4),y=c(0.25,0.5)),
                                    aspectmode='data'),
                      scene8 = list(domain=list(x=c(0.3,0.7),y=c(0.25,0.5)),
                                    aspectmode='data'),
                      scene9 = list(domain=list(x=c(0.6,1),y=c(0.25,0.5)),
                                    aspectmode='data'),
                      scene10 = list(domain=list(x=c(0.3,0.7),y=c(0.,0.25)),
                                    aspectmode='data'),
                      scene11 = list(domain=list(x=c(0.6,1),y=c(0.,0.25)),
                                    aspectmode='data'))

fig

Now we can run a few Herring TS simulations. We wil use 1 random length for all the shapes and simulate for 38 kHz at a rotation angle of 90°:

set.seed(78965);L=rnorm(1,14,2.3) #random lengths drawn from a normal distribution with a mean of 14 cm and an sd of 2.3 cm

ts.all=data.frame()
for(i in 1:11){
  message(Sys.time(),": Processing Shape ", i)
  ts=krm(frequency = c(38)*1e3,
      c.w = 1490,
      rho.w = 1030,
      theta = 90,
      cs = c(1570, 345),
      rhos = c(1070, 1.24),
      shape = readRDS(paste0("HER_",sprintf("%02d", i),"_shape.rds")),
      modes = c("fluid", "soft"),
      L = L/100,
      fb = 1)
  ts$shape=i
  ts.all=rbind(ts.all,ts)
}
## 2023-03-06 20:55:11: Processing Shape 1
## 2023-03-06 20:55:11:2 Body parts detected with 2 densities and 2 internal soundspeeds
## 2023-03-06 20:55:11: Processing Shape 2
## 2023-03-06 20:55:11:2 Body parts detected with 2 densities and 2 internal soundspeeds
## 2023-03-06 20:55:11: Processing Shape 3
## 2023-03-06 20:55:11:2 Body parts detected with 2 densities and 2 internal soundspeeds
## 2023-03-06 20:55:11: Processing Shape 4
## 2023-03-06 20:55:11:2 Body parts detected with 2 densities and 2 internal soundspeeds
## 2023-03-06 20:55:11: Processing Shape 5
## 2023-03-06 20:55:11:2 Body parts detected with 2 densities and 2 internal soundspeeds
## 2023-03-06 20:55:11: Processing Shape 6
## 2023-03-06 20:55:11:2 Body parts detected with 2 densities and 2 internal soundspeeds
## 2023-03-06 20:55:11: Processing Shape 7
## 2023-03-06 20:55:11:2 Body parts detected with 2 densities and 2 internal soundspeeds
## 2023-03-06 20:55:11: Processing Shape 8
## 2023-03-06 20:55:11:2 Body parts detected with 2 densities and 2 internal soundspeeds
## 2023-03-06 20:55:11: Processing Shape 9
## 2023-03-06 20:55:11:2 Body parts detected with 2 densities and 2 internal soundspeeds
## 2023-03-06 20:55:11: Processing Shape 10
## 2023-03-06 20:55:11:2 Body parts detected with 2 densities and 2 internal soundspeeds
## 2023-03-06 20:55:11: Processing Shape 11
## 2023-03-06 20:55:11:2 Body parts detected with 2 densities and 2 internal soundspeeds
kable(ts.all)
frequency c.w rho.w theta L scale f1 f2 sigma TS shape
38000 1490 1030 90 0.1336088 0.1336088 0.0010120-0.0000497i -0.0000083+0.0068226i 0.0068468 -43.29020 1
38000 1490 1030 90 0.1336088 0.1336088 0.0014130-0.0011905i -0.0014107+0.0079948i 0.0068043 -43.34435 2
38000 1490 1030 90 0.1336088 0.1336088 0.0011609-0.0006970i -0.0021915+0.0059686i 0.0053714 -45.39824 3
38000 1490 1030 90 0.1336088 0.1336088 0.0002514-0.0007578i 0.0018561+0.0106006i 0.0100659 -39.94298 4
38000 1490 1030 90 0.1336088 0.1336088 -0.0000246-0.0006323i -0.0020224+0.0076785i 0.0073375 -42.68904 5
38000 1490 1030 90 0.1336088 0.1336088 0.0005926-0.0006744i -0.0006447+0.0089434i 0.0082691 -41.65086 6
38000 1490 1030 90 0.1336088 0.1336088 0.0007244-0.0007589i 0.0010252+0.0108122i 0.0102045 -39.82420 7
38000 1490 1030 90 0.1336088 0.1336088 0.0008308-0.0005894i -0.0012616+0.0076586i 0.0070823 -42.99653 8
38000 1490 1030 90 0.1336088 0.1336088 0.0002398-0.0009108i 0.0001149+0.0109894i 0.0100849 -39.92653 9
38000 1490 1030 90 0.1336088 0.1336088 0.0009510-0.0014875i -0.0012592+0.0094089i 0.0079274 -42.01739 10
38000 1490 1030 90 0.1336088 0.1336088 0.0013105-0.0006044i -0.0039944+0.0061406i 0.0061525 -44.21902 11